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To  explore  potential  for  the  power  density  enhancement  of  solid  oxide  fuel  cells  by  controlling  the 
cathode-electrolyte  interface  in  mesoscale,  two-dimensional  numerical  simulations  were  conducted. 
In  the  simulation,  a  level  set-based  topology  optimization  technique  was  successfully  coupled  with  the 
SOFC  simulation  based  on  a  microscale  model  and  was  applied  for  the  local  optimization  of  the  interface 
shape.  The  numerical  results  showed  that  the  optimized  shape  of  the  cathode-electrolyte  interface  varied 
depending  on  the  simulation  conditions  and  that  the  cell  performance  could  be  improved  by  applying 
non-flat  design  to  the  cathode-electrolyte  interface  for  the  same  amount  of  cathode/electrolyte  materials. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  devices  which  directly  convert  the  chemical  energy 
of  reactants  into  electricity  through  electrochemical  reactions. 
Among  various  fuel  cell  types,  the  solid  oxide  fuel  cell  (SOFC) 
has  recently  been  attracting  particular  attention  for  its  high  effi¬ 
ciency  and  fuel  flexibility  [1].  One  of  the  unique  features  of  SOFCs 
is  that  all  components  in  an  SOFC  system  are  solid,  including  the 
ceramic  electrolyte,  which  permits  exceptional  geometrical  flexi¬ 
bility  of  the  cell  design.  By  controlling  the  cell’s  macroscopic  shape, 
we  can  achieve  designs  that  offer  select  advantages.  For  exam¬ 
ple,  a  tubular  type  of  cell  has  good  resistance  to  thermal  stress 
while  planar,  or  flattened-tubular  types,  can  be  stacked  compactly 
to  provide  high  power  density  [2,3].  A  micro-tubular  type  of  cell 
is  mechanically  strong  and  can  achieve  the  high  power  densities 
suitable  for  small  applications  [4,5].  Geometry  control  is  a  chal¬ 
lenging  topic  when  designing  the  porous  electrodes  of  a  SOFC.  It  is 
widely  recognized  that  the  electrode  microstructure  has  a  signifi¬ 
cant  impact  on  cell  performance  as  well  as  cell  durability  [6-10]. 
Therefore,  quantification  of  the  porous  microstructure  and  its  cor¬ 
relation  to  the  cell  performance  are  the  recent  interests  [11,12]. 
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These  porous  electrodes  must  have  appropriate  microstructures 
that  increase  the  density  of  reaction  sites,  reduce  resistance  to 
gas  diffusion  and  establish  efficient  networks  of  both  ion  and 
electron  conductors.  Unlike  the  direct  control  of  geometry  in  the 
macroscale,  the  electrode  microstructure  is  usually  controlled  indi¬ 
rectly  and  statistically.  Composite  electrode  (e.g.,  LSM-YSZ)  is  an 
example  of  microstructure  control  commonly  used  in  SOFC  [13]. 
Following  recent  understanding  that  electrochemical  reaction  pro¬ 
ceeds  most  prominently  near  the  electrode-electrolyte  interface, 
the  microstructure  control  within  this  effective  thickness  of  the 
electrode  is  gathering  attention.  Reports  on  the  surface  roughness 
of  the  interface  can  be  found  in  literature  [14,15].  Use  of  nanopar¬ 
ticles  near  the  interface  to  increase  the  number  of  reaction  sites  is 
also  an  example  of  a  microstructure  modification  [16-18].  How¬ 
ever,  this  approach  currently  encounters  durability  problems  in 
long-term  use  [9,10,19]. 

Konno  et  al.  [20]  recently  investigated  the  possibility  of  the 
electrode-electrolyte  interface  shape  modification  to  increase  the 
interface  area  and  consequently  the  active  reaction  region  in 
the  electrode.  The  order  of  the  structure  scale  in  the  study  was 
about  10-100  pum  so  that  the  characteristic  length  of  the  structure 
became  larger  than  the  effective  thickness  of  the  electrode.  The 
effects  of  this  “mesoscale”  structure  on  Ni-YSZ  anode  performance 
were  examined  numerically  and  experimentally  and  positive  con¬ 
clusions  were  obtained  on  this  concept.  It  is  worth  noting  that  a 
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Nomenclature 

Dy  binary  diffusivity  (m2  s-1 ) 

DiK  Knudsen  diffusivity  (m2  s-1 ) 

d  particle  diameter  (m) 

E0  EMF  at  the  standard  state  (V) 

F  Faraday  constant  (C  mol-1 ) 

Itpb  charge  transfer  current  (A  m-3 ) 

i  current  density  (A  nrr2 ) 

itpb  charge  transfer  current  per  TPB  length  (A  nrr 1 ) 

I<  permeability 

ltpb  TPB  density  (nrr2) 

P  percolation  probability 

p  partial  pressure  (Pa) 

R  universal  gas  constant  (J  mol-1  K-1 ) 

T  temperature  (K) 

V  volume  fraction 

Vmax  upper  limit  of  volume  constraint  (m3 ) 

Greek  letters 

s  porosity  of  porous  electrode 

0  potential  (V) 

0o  Nernst  potential  (electromotive  force)  (V) 

0fc  cell  terminal  voltage  (V) 

y  coefficient  of  complexity 

rj  overpotential  (V) 

a  conductivity  (S  nr1 ) 

r  tortuosity  factor 

Subscripts 
A  anode 

act  activation 

ave  average 

C  cathode 

E  electrolyte 

el  electronic  conducting  material 

io  ionic  conducting  material 

tpb  three-phase  boundary 

oo  bulk  gas  phase 


geometry  control  of  this  scale  has  some  advantages  over  macro¬ 
scopic  or  microscopic  scales.  Because  the  structure  is  large  enough 
compared  to  the  microstructure  of  the  porous  electrode,  it  is  free 
from  the  durability  problem  caused  by  the  morphology  change.  Pre¬ 
cision  fabrication  of  the  designed  interface  shape  is  possible  as  well. 
On  the  other  hand,  the  structure  is  small  enough  compared  to  the 
macroscopic  cell  shape  that  it  can  be  applied  in  a  variety  of  cell 
types  (e.g.,  tubular  and  planar). 

An  effective  design  of  the  electrode-electrolyte  interface  is  not 
clear  yet  because  the  effectiveness  of  the  interface  shape  control  in 
mesoscale  depends  on  various  resistances  in  the  cell  and  the  oper¬ 
ation  condition.  For  such  design  problem,  there  is  a  possibility  that 
an  optimization  calculation  technique  can  be  efficiently  used.  In 
this  paper,  to  explore  the  potential  for  power  density  enhancement 
through  mesoscale  interface  control,  we  perform  two-dimensional 
numerical  simulations  of  a  SOFC  where  a  level  set-based  topology 
optimization  technique  [21  ]  is  applied  for  the  local  optimization  of 
the  interface  shape. 

2.  Concept  of  shape  modification  of  the 
electrode-electrolyte  interface 

The  electrochemical  reactions  proceed  at  the  three-phase 
boundary  (TPB)  where  the  ion  conductor,  electron  conductor  and 


gas  phase  meet.  TPBs  are  distributed  in  the  porous  electrodes  but  it 
is  reported  that  they  are  not  equally  active,  that  the  electrochemical 
reactions  occur  more  prominently  near  the  electrode-electrolyte 
interface  [22,23].  Thus  the  most  effective  thickness  of  electrodes  is 
an  open  question,  but  in  this  section,  for  simplicity,  we  assume  that 
the  electrochemical  reactions  take  place  only  at  the  interface. 

Fig.  1  schematically  shows  the  concept  of  an  interface  modifi¬ 
cation,  taking  the  cathode  side  as  an  example.  If  the  interface  has 
a  wavy  shape  as  shown  in  Fig.  1(b),  the  interface  area  is  increased 
compared  to  the  flat  case  (Fig.  1(a)),  which  implies  that  the  elec- 
trochemically  active  region  is  enlarged  by  the  wavy  shape.  At  the 
position  marked  “A”  in  Fig.  1(b),  the  electrolyte  thickness  is  reduced 
while  the  cathode  thickness  is  increased.  This  means  that  the  poten¬ 
tial  loss  associated  with  ion  transfer  in  the  electrolyte  is  reduced, 
but  the  concentration  overpotential  is  increased.  At  position  “B”,  an 
opposite  situation  is  obtained.  The  overall  performance  of  the  cell 
is  determined  by  the  sum  of  all  these  effects. 

We  further  simplify  the  problem  by  neglecting  the  concentra¬ 
tion  overpotential  because  its  effect  is  usually  minor  in  a  practical 
cell  operating  under  normal  conditions.  Then  the  performance  of  a 
flat  cell  (Fig.  1(a))  can  be  described  by  the  balance  between  the  ion 
transfer  resistance  in  the  electrolyte  and  the  reaction  resistance 
at  the  interface.  Konno  et  al.  [20]  introduced  a  non-dimensional 
number,  BiFC,  as  follows: 

gj  _  ASR electrolyte  _  1  / ASRreaction  _  1  / ASRreactjon  ^ 

ASRreaction  1  /  ASReiectrolyte  &  electrolyte  /^-electrolyte 

ASR,  <y electrolyte* and  leiectroiyte are  the  area  specific  resistance,  ion  con¬ 
ductivity  of  the  electrolyte,  and  electrolyte  thickness,  respectively. 
ASRreaction  can  be  expressed  in  terms  of  the  activation  overpotential 
and  the  current  density. 

ASRreaction  =  ''f  (2) 

When  BiFc  is  larger  than  unity,  ion  transport  is  the  dominant 
phenomenon  and  expansion  of  the  interface  area  becomes  ineffec¬ 
tive.  Expansion  of  the  interface  area  can  be  effective  when  BiFC  is 
less  than  unity,  where  the  reaction  resistance  is  the  bottleneck.  The 
value  of  BiFc  is  affected  by  the  choice  of  electrolyte  materials,  elec¬ 
trolyte  thickness,  temperature,  current  density  and  microstructure 
of  the  electrode.  We  assume  a  cell  made  of  Yttria  Stabilized  Zirco- 
nia  (YSZ)  electrolyte  and  Lanthanum  Strontium  Manganate  (LSM) 
cathode  and  estimate  the  value  of  BiFC.  creiectroiyte  and  pact  are  esti¬ 
mated  following  Bessette  et  al.  [24]  and  Costamagna  and  Flonegger 
[25],  respectively.  BiFC  values  for  different  electrolyte  thicknesses 
and  temperatures  are  shown  in  Fig.  2.  The  current  density  is  fixed 
at  300  mA cm-2.  BiFC  can  assume  values  less  than  1,  especially 
when  the  electrolyte  is  thin.  Considering  that  the  electrolyte  thick¬ 
nesses  of  today’s  electrode-supported  cells  are  typically  in  the 
range  of  10-20  p,m,  Fig.  2  indicates  the  possibility  of  power  density 
enhancement  by  shape  modification  of  the  electrode-electrolyte 
interface. 

The  above  discussion  by  BiFC  is  based  on  the  assumption  that 
the  electrochemical  reactions  proceed  at  the  electrode-electrolyte 
interface.  The  effect  of  concentration  overpotential  is  also 
neglected.  In  reality,  the  cell  performance  is  affected  by  the  ionic 
resistivity  in  electrolyte,  ionic  and  electronic  resistivity  in  elec¬ 
trode,  electrode  reaction  resistivity  and  gas  diffusion.  To  consider 
these  effects,  we  perform  numerical  simulation  and  the  results  are 
reported  in  the  following  sections.  Although  interface  modifica¬ 
tion  is  expected  to  be  particularly  effective  when  the  electrolyte 
is  thin  as  shown  in  Fig.  2,  electrolyte-supported  cells  with  thick 
electrolytes  are  used  in  this  study.  BiFC  of  the  cell  assumed  in 
the  following  simulation  is  about  4,  so  the  ion  transfer  resistance 
in  the  electrolyte  is  the  major  resistance  of  the  cell.  By  using  an 
electrolyte-supported  cell,  we  can  limit  the  shape  modification 
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(a)  Flat  interface 


Electrolyte 


Image  of  '  LSM 
Microstructure 


(b)  Wavy  interface 


Fig.  1.  Image  of  electrolyte-electrode  interface  transformation. 


to  only  one  side  of  the  cell  (the  cathode)  without  affecting  the 
shape  of  the  other  side  (the  anode).  This  is  ideal  for  discussion 
in  this  study  because  we  can  avoid  simultaneous  modification  of 
the  two  interfaces,  namely,  the  cathode-electrolyte  interface  and 
the  anode-electrolyte  interface.  It  should  be  noted  that  our  goal  in 
this  study  is  not  to  develop  high  power  density  cells,  but  to  obtain 
fundamental  information  on  how  interface  modification  affects  the 
cell’s  power  density,  and  to  clarify  if  this  is  a  promising  technique 
for  developing  more  advanced  SOFCs. 

3.  Numerical  procedure 

Fig.  3  shows  a  schematic  view  of  the  cell  considered  in  this 
study.  It  is  an  electrolyte-supported  cell,  and  the  relatively  thick 
electrolyte  is  sandwiched  between  anode  and  cathode  electrodes. 
The  electrolyte,  anode  and  cathode  materials  are  assumed  to  be 
YSZ,  Ni-YSZ  cermet  and  LSM-YSZ  composite,  respectively.  In  this 
study,  an  optimization  calculation  is  implemented  into  the  SOFC 
power  generation  simulation.  The  computational  domain  for  the 
SOFC  power  generation  covers  all  anode,  electrolyte  and  cathode 
layers  while  the  computational  domain  for  the  optimization  cal¬ 
culation  is  limited  to  a  half  cell  region  including  only  the  cathode 
and  upper  half  of  the  electrolyte,  as  shown  in  Fig.  3.  The  width  of 
the  computational  domain  is  set  to  50  pm  to  appropriately  limit  the 
computational  load.  The  following  assumptions  were  made  for  sim¬ 
plification.  The  problem  was  assumed  to  be  two-dimensional  and 
steady.  An  isothermal  condition  was  applied,  assuming  a  constant 


temperature  of  800  °C.  The  microstructure  of  the  porous  electrodes 
was  modeled  by  the  random  packing  of  spheres  with  a  constant 
and  uniform  particle  diameter  for  both  ionic  and  electronic  con¬ 
ductors  in  the  electrodes.  Based  on  these  assumptions,  a  numerical 
calculation  program  was  developed.  The  next  section  describes  an 
overall  picture  of  the  numerical  model  for  the  SOFC  power  gener¬ 
ation  prediction,  followed  by  a  section  describing  the  optimization 
method  applied  to  the  cathode  side  of  the  cell  in  this  study. 

3.1.  Simulation  of  SOFC  power  generation 

3.1.1.  Electrode  microscale  model 

We  adopted  the  microscale  model  proposed  by  Nam  and  Jeon 
[26]  and  extended  it  to  a  2D  simulation.  This  model  calculates 
the  charge  transfer  current  at  the  three-phase  boundary  from  the 
volume  and  area  specific  lengths  of  the  three-phase  boundary 
evaluated  by  a  model  based  on  random  packing  of  spheres,  and 
empirical  correlations  for  the  charge  transfer  current  density  based 
on  experiments  with  patterned  mesh  electrodes.  The  anode  and 
cathode  charge  transfer  current  density  for  a  unit  length  of  the 
three-phase  boundary  proposed  in  Ref.  [26]  are  derived  from  the 
experiments  of  Bieberle  et  al.  [27]  and  Radhakrishnan  et  al.  [28], 
respectively,  and  expressed  as 

('kpb'>A  ~  1.645 Ph°  "pS°o7  exp(10,  212 /T)  x  10"4’'"  ’  ^ 


Temperature  [°C] 


l/202+2e-— »Cr 


f  Computational  domain  for 
optimization  calculation 


Cathode 

T 

j  50pm 

i 

Electrolyte 

| 

50|im 

200pm 

i. 

Anode 

\  50  pm 

H2+0“’— KH2C)+2e- 


*-  Computational  domain  for 
SOFC  calculation 


Fig.  2.  BiFC  estimation  for  YSZ  electrolyte  and  LSM  cathode  system  at  an  average  Fig.  3.  Computational  domains  for  the  SOFC  calculation  and  for  the  optimization 
current  density  of  300  mA  cm-2 .  calculation. 
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Table  1 

Geometric  parameters  of  electrode  microscale  structure. 


Value 

0.3 

ta,  rc 

3 

(Vel)A,(Vei)c 

0.5 

C Vi0)A,(Vi0)c 

0.5 

( nei)A ,  ( nd)c 

0.5 

(ni0)A,  (ni0)c 

0.5 

(dei) a,  (del)c 

1  |jim 

( dio )/l ,  (dio)c 

1  jxm 

ec 

30° 

RT  2  sinh[2Fr/gct/RT)  (4) 

tpb  c  4F  o.OOl 36p-°'25exp(l 7,401  /T)’ 

where  rjact  is  the  activation  overpotential  and  p,-  is  the  partial  pres¬ 
sure  of  species  i.  The  charge  transfer  current  is  calculated  as 


hpb  —  hpbhpbi 


(5) 


where  ltpb  is  the  density  of  the  three-phase  boundary  per  unit  vol¬ 
ume.  It  is  evaluated  by  the  random  packing  model  using  values  for 
porosity,  e\  tortuosity  factor,  r;  volume  ratio  of  ion  and  electron 
conductors,  VehVio;  number  density  of  ion  and  electron  conduc¬ 
tors,  nei,nio ;  particle  diameter  of  ion  and  electron  conductors,  de/,d10 ; 
contact  angle  of  the  particles,  6C\  as  shown  in  Table  1. 

The  activation  overpotential  is  calculated  as 


RT 

{if act) a  —  (0e/  -  0io)  +  ^ 


(  Ph2  Ph2Q,oo 
\PH2,oo  Ph2o 


(6) 


tooctfc  =  (&,-**) +5  lnf^)  (7) 

\Po2,oo  J 

The  last  term  of  the  right-hand  side  of  Eqs.  (6)  and  (7)  corresponds 
to  the  concentration  overpotential.  Equations  (3),  (4),  (6)  and  (7) 
determine  the  electrochemical  reaction  rates  for  a  given  condition 
including  the  ionic  and  electronic  potentials  and  gas  concentra¬ 
tions. 


3 A. 2.  Charge  transfer  and  gas  diffusion  model 

The  ionic  and  electronic  conductivities  (S  cm-1 )  of  the  electrode 
and  electrolyte  bulk  materials  are  obtained  following  expressions 
reported  in  literature  [24,26,29,30]. 


<z=(0.00294exp(10f°))-1 

(8) 

o  4.2  xlO5  /-1200\ 

ctlsm  =  T  exp(  T  ) 

0) 

cr°;  =  3.27  x  104  -  10.653T 

(10) 

The  effective  ionic  or  electronic  conductivities  in  the  porous  elec¬ 
trode  are  expressed  as 

ai=a°[(\-e)ViPi]h5, 

(11) 

where  P*  is  the  percolation  probability  calculated  through  the  ran¬ 
dom  packing  model.  The  effective  conductivities  are  summarized 
in  Table  2.  The  effective  electronic  conductivity  is  significantly  large 
compared  to  the  effective  ionic  conductivity.  Therefore,  we  assume 
a  uniform  electronic  potential  in  both  electrodes  and  consider  only 


Table  2 

Effective  conductivity  of  electrodes. 


Anode 

Cathode 

oi0  (Scrrr1) 

4.1  x  10-3 

4.1  x  10“3 

crel  (Scrrr1) 

3.6  xlO3 

2.4  xlO2 

the  transport  of  ions  and  mass  in  the  following  computational 
model. 

Fig.  4  schematically  shows  the  equivalent  circuit  model  in  one¬ 
dimensional  geometry.  Both  ion  transport  and  gas  diffusion  are 
considered  in  the  electrodes  while  only  ion  transport  takes  place 
in  the  electrolyte.  The  governing  equation  for  ion  transport  is 
expressed  as 


B  = 


-1  Anode 
0  Electrolyte 
+1  Cathode 


(12) 


The  term  in  the  right-hand  side  is  associated  with  the  charge  trans¬ 
fer  between  the  ion  conductive  phase  and  the  electron  conductive 
phase.  It  is  evaluated  using  the  microscale  model  shown  in  the 
previous  section. 

For  the  mass  transport  in  the  anode,  the  dusty-gas  model  is 
adopted  [31-33].  This  model  includes  the  diffusion  flow  and  the 
viscous  flow  and  is  derived  from  Maxwell-Stefan  equations  consid¬ 
ering  the  effect  of  Knudsen  diffusion  and  Darcy’s  law.  The  governing 
equations  can  be  expressed  as 


d  f  /<h2  9Ph2  \  d  f  fh2  dpn2  \  d  /  ^<dPh2  d pt  \ 
dx  \W~dx~ J  +  dy  \W~dy~ J  +  dx  \  ^T~dx  J 

d  f  /<dPh2  dpt  A  _  hpb 
+  dy  RT  dy  J  ~  2F’ 


(13) 


d  (  /<h2o  9ph2o  \  d  f  /<h2o  9pn2o  \  d  (  1<dPh2o  dpt  \ 

dx  yW  dx  J  +  dy  [fW  dy  J  +  dx  \  RT  Ik  J 

d  f  1<dPh2o  dpt\  _  hpb 

9y  V  RT  &  J  ~  2F 


l<  H2  = 


D 


0 


D 


0 


H2-H2O^H2,K 


D\ 


0 

H2-H20 


+  D 


0  ’ 
m,I< 


1<h2o  = 


D 


0 


0 


h2-h2o^h2o  ,k 


D 


0 

H2-H20 


+  D\ 


0 

m,K 
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Fig.  5.  Flowchart  of  the  optimization  calculation  coupled  with  the  SOFC  simulation. 
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governing  equation  is  expressed  as 
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The  effective  diffusivities  for  binary  diffusion  and  Knudsen  diffu¬ 
sion  are  affected  by  properties  of  the  porous  electrodes  and  are 
evaluated  as 
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r  is  the  tortuosity  factor  of  porous  electrodes.  D2_j  and  DlK  are  the 
binary  and  Knudsen  diffusivity,  respectively. 


where  pt  is  total  pressure,  x2  is  the  molar  fraction,  and  /z  is  the 
viscosity  of  the  gas  mixture.  I<  is  the  permeability  of  the  porous 
electrode,  which  is  evaluated  using  the  Carman-Kozeny  relation 
[34]  as 


5(1  -  e)2(6/d)2 

where  d  is  the  particle  diameter  of  the  ionic/electronic  conductors. 
The  mass  transport  in  the  cathode  is  evaluated  as  a  self-diffusion 
of  oxygen  assuming  the  molar  flux  of  nitrogen  is  zero  [35].  The 


3.2.  Level  set-based  topology  optimization 
3.2 A.  Formulation 

The  optimization  method  employed  in  this  study  is  briefly  intro¬ 
duced  in  this  section.  Details  of  the  method  can  be  found  in  Ref.  [21  ]. 
D  is  a  fixed  design  domain  consisting  of  the  electrolyte  and  cathode 
and  the  corresponding  area  is  marked  in  Fig.  2.  We  define  the  level 
set  function,  t/f(x),  as 
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Fig.  6.  Evolution  of  the  cathode-electrolyte  interface  during  optimization  calculations. 
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Fig.  8.  Influence  of  the  initial  configuration  upon  optimization  results. 


Positive  values  of  ifr(x)  correspond  to  the  cathode  area,  negative 
values  to  the  electrolyte  and  the  interface  between  the  electrolyte 
and  the  cathode  is  expressed  by  lines  where  i/r(x)  =  0.  The  objec¬ 
tive  functional,  F,  and  the  constraint  functional  with  respect  to  the 
volume,  G,  are  defined  as 


rm))=  1 

f  f(x)d£2 

(19) 

JQ 

cm))=  . 

1  dQ  —  Vmax  <  0 

(20) 

JQ 

Vniax  is  the  upper  limit  of  the  volume  constraint.  We  replace  this 
optimization  problem  with  an  optimization  problem  to  minimize 
the  energy  functional,  which  is  the  sum  of  the  fictitious  interface 
energy  and  the  objective  functional,  formulated  as  follows: 

FRm),f)=  [  f(x)dC2+  [  ly\Vxfr\2dC2  (21) 

Jn  Jd  z 

G(n(1r))=  [  d£2  —  1/max  <  0  (22) 

JC2 

where  FR  is  a  regularized  objective  functional  and  y  is  a  regular¬ 
ization  parameter  representing  the  ratio  of  the  fictitious  interface 


energy  and  the  objective  functional.  The  regularization  parame¬ 
ter  is  possible  to  adjust  the  degree  of  complexity  of  the  optimized 
structures  and  sometimes  called  as  the  “complexity  parameter”. 

The  optimization  problem  represented  by  Eqs.  (21)  and  (22)  is 
replaced  with  an  optimization  problem  without  constraints  using 
Lagrange’s  method  of  undetermined  multipliers.  That  is,  when  the 
Lagrangian  is  defined  as  FR  and  the  Lagrange  multiplier  of  the  vol¬ 
ume  constraint  is  k,  the  optimization  problem  represented  by  Eqs. 
(21 )  and  (22)  is  replaced  with 
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We  also  assume  that  level  set  function  changes  are  proportional 
to  the  gradient  of  Lagrangian  FR,  based  on  the  concept  of  the  phase 
field  method,  as  shown  in  the  following, 


dxjf 

~dt 


-KMjf 


(24) 


where  K(ijf)  >  0  is  a  constant  of  proportion.  Substituting  Eq.  (23)  into 
Eq.  (24),  we  have 


^  =-K(ir)(DTF-YV2ir) 


(25) 


where  DTF  is  a  topological  derivatives  with  respect  to  the 
Lagrangian  FR.  The  Lagrange  multiplier  is  expressed  as 


„  _  f^KW(DTF  +  yV2x/f)dF2 
JQm)dF2 


(26) 


To  represent  the  level  set  function  independent  of  the  exterior  of 
the  fixed  design  domain  D,  we  assume  that  the  boundary  condition 
of  the  level  set  function  obeys  a  Dirichlet  boundary  condition  on 
the  non-design  boundary,  and  a  Neumann  boundary  condition  on 
the  other  boundaries.  The  system  of  the  time  evolutional  equation 
can  then  be  represented  as 
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Fig.  9.  Influence  of  O2  diffusivity  upon  optimization  results,  y  =  0.02. 
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Fig.  10.  Ion  potential,  activation  overpotential,  concentration  overpotential  and  three-phase  boundary  charge  transfer  current  for  optimal  configuration,  kdkcstd  =  1 .0,  y  =  0.02. 


3.2.2.  Application  to  SOFC  simulation 

The  design  domain  includes  the  cathode  and  upper  half  of  the 
electrolyte,  as  shown  in  Fig.  3.  As  mentioned  above,  a  Dirichlet 
boundary  condition  is  set  at  the  cathode  surface  while  a  Neumann 
boundary  condition  is  applied  to  the  other  boundaries.  The  appar¬ 
ent  average  current  density  is  expressed  as 


In  hpbhpbdfi 

$ cathode 


(28) 


S cathode  is  the  cathode  surface  area  within  the  computational 
domain.  We  set  the  objective  functional  F(£2{\jr))  as  below  and  con¬ 
sider  an  optimization  problem  to  maximize  the  total  current  under 
a  given  terminal  voltage. 


0  —  hpbhpb  (29) 

Fm 0)  =  -  [  itpbltpbdV  (30) 
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The  Lagrangian  is  expressed  as 

F=  f  (f(x)  +  h)d£2  =  f  (-itpbltpb  +  ^)d£2 

J  £2  J  Q 

(31) 

Ionic  conductivity  and  the  gas  diffusion  coefficient  are  given  as 

<*io  =  +  Oe(1  -  H(f)) 

(32) 

ko2  =  kcHW) 

(33) 

H(xjA  _  j  1  at  cathode 

1  v ;  0  at  electrolyte 

(34) 

where  crc  and  op  are  the  effective  ionic  conductivity  in  the  cath¬ 
ode  and  the  ionic  conductivity  of  the  electrolyte,  respectively.  The 
maximum  volume  of  the  cathode  is  set  so  that  the  average  cath¬ 
ode  thickness  does  not  exceed  50  p,m.  This  is  equivalent  to  ensure 
the  minimum  average  electrolyte  thickness  of  200  |jim.  To  avoid  an 
artificial  formation  of  electrolyte  inside  the  cathode,  the  following 
condition  is  applied  inside  the  cathode. 

/(*)  =  0  if  /(*)<  0  (35) 

The  flowchart  of  the  simulation  is  shown  in  Fig.  5. 
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Fig.  11.  Influence  of  the  complexity  parameter  upon  optimization  results  (/<c//<c_std  =  0.4). 
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3.3.  Computational  conditions 


4.2.  Optimization  under  different  conditions 


The  cell  is  assumed  to  operate  under  atmospheric  pressure  at 
a  constant  and  uniform  operation  temperature  of  800  °C.  The  fuel 
supplied  to  the  anode  surface  is  humidified  hydrogen,  the  gas  com¬ 
position  of  which  is  H2  :H20  =  0.99:0.01 .  On  the  cathode  surface,  dry 
oxygen  diluted  with  nitrogen  is  supplied.  The  cathode  gas  com¬ 
position  is  02:N2  =  0.05:0.95.  The  anode  and  cathode  electronic 
potentials  are  given  as  the  boundary  conditions  for  calculations 
of  the  equivalent  circuit.  We  applied  the  Nernst  potential  to  the 
anode  while  the  cell  terminal  voltage  was  applied  to  the  cathode. 
The  Nernst  potential  is  calculated  as 


RT  . 

'~Eo+ 2 F  n 


(  PH2,oc  /  Po2,oc  \°'5\ 
Vph2o,ooUo,i3o;  ) 


(36) 


where  E0  is  the  electromotive  force  at  the  standard  state.  Symmetric 
boundary  conditions  are  applied  to  the  vertical  boundaries  in  the 
x-direction  in  Fig.  3. 


4.  Simulation  results  and  discussion 

4.1.  General  features  of  optimization  calculation 

4.1  A.  Example  of  optimizing  process 

Fig.  6  shows  a  typical  example  of  the  shape  evolution  dur¬ 
ing  optimization  calculations.  The  thick  gray  region  represents 
the  cathode  area.  Changes  in  the  average  current  density  and 
Lagrangian  during  these  particular  calculations  are  shown  in  Fig.  7. 
The  step  number  in  the  figures  indicates  the  number  of  times  the 
level  set  function  has  been  updated.  Step  0  is  the  initial  state.  Right 
after  the  calculation  starts,  the  dividing  line  between  the  cath¬ 
ode  and  the  electrolyte  gradually  moves  as  shown  in  Fig.  6.  The 
volume  of  the  cathode  increases  until  the  step  number  reaches 
approximately  100  when  a  sudden  change  in  the  Lagrangian  is 
observed,  as  shown  in  Fig.  7.  This  is  when  the  cathode  reaches 
the  upper  limit  of  the  volume  constraint.  The  dividing  line  con¬ 
tinues  to  evolve  after  step  number  100  but  the  volume  of  the 
cathode  (or  electrolyte)  is  kept  constant.  The  shape  of  the  inter¬ 
face  is  almost  stable  after  step  number  1200,  corresponding  well 
with  the  almost  constant  average  current  density  after  step  1200 
observed  in  Fig.  7.  The  average  current  density  is  increased  from 
its  initial  value  of  678  mA cm-2  to  the  final  value,  726mA cm-2, 
during  the  calculation  process,  showing  that  the  level  set-based 
topology  optimization  method  is  successfully  implemented  in  this 
2D  SOFC  calculation.  To  the  best  of  the  authors’  knowledge,  this 
is  the  first  report  of  an  optimization  calculation  applied  in  SOFC 
simulations. 


4.1.2.  Effect  of  initial  configuration 

To  clarify  if  the  results  show  any  dependency  on  the  initial 
configuration,  optimization  calculations  are  conducted  for  the  dif¬ 
ferent  initial  configurations  shown  in  Fig.  8.  The  only  difference 
among  these  calculations  is  the  initial  configuration;  all  other  con¬ 
ditions  are  identical.  The  final  results  obtained  by  the  optimization 
calculations  match  each  other  very  well,  as  can  be  confirmed  in 
Fig.  8.  Note  that  the  final  result  shown  in  Fig.  8(c)  can  be  con¬ 
sidered  practically  the  same  as  that  of  Fig.  8(a)  and  (b)  because 
a  symmetric  boundary  condition  is  applied  in  the  x-direction  for 
these  calculations.  The  average  current  density  at  the  final  states 
are  725.60  mA  cm-2  in  all  cases.  It  demonstrates  that  the  present 
optimization  calculation  method  is  independent  of  the  initial  con¬ 
figuration. 


The  effects  of  cathode-electrolyte  interface  modification  on  the 
cell’s  power  generation  performance  are  discussed  in  this  section, 
based  on  the  numerical  results  obtained  by  the  optimization  cal¬ 
culations.  Calculations  were  conducted  for  different  values  of  the 
oxygen  diffusivity  and  complexity  parameter. 

4.2. 1 .  Effect  of  oxygen  diffusivi ty 

The  value  of  the  oxygen  diffusivity,  kc,  which  appears  in  Eq. 
(16),  is  artificially  changed  to  investigate  how  the  gas  diffusion 
affects  the  optimized  geometry.  The  oxygen  diffusivity  estimated 
through  the  model  explained  in  Section  3.1  is  set  as  the  standard 
value,  expressed  here  as  kc_std-  Calculations  are  carried  out  for  dif¬ 
ferent  values  of  the  oxygen  diffusivity,  with  kclK-std  =  1  -0,  0.4  and 
0.1.  In  a  real  cathode,  a  change  in  kc  should  be  associated  with 
some  change  of  the  porous  microstructure  and  therefore  other  geo¬ 
metric  parameters  are  expected  to  change  simultaneously.  In  this 
study,  however,  such  effects  are  neglected  and  only  the  value  of  kc 
is  altered. 

Fig.  9  shows  the  shapes  of  the  optimized  cathode-electrolyte 
interface  for  the  three  cases  calculated.  It  is  obvious  that  the 
value  of  the  oxygen  diffusivity  affects  the  final  results.  A  deep 
and  large-scale  groove  is  observed  when  /<c/kc_std  =  1  -0.  The  grooved 
part  becomes  shallower  for  smaller  values  of  kclkcj;td  and  when 
kclkc_std  =  0.1,  the  interface  is  almost  flat. 

To  evaluate  power  generation  enhancement,  we  conducted  a 
calculation  for  a  reference  case  in  which  both  the  cathode  and 
electrolyte  shapes  are  set  as  flat,  with  constant  and  uniform  thick¬ 
nesses  of  50  |jim  and  200  |jim,  respectively.  The  value  of  kclkCJStd  is 
unity.  The  other  conditions  are  the  same  as  for  the  optimization 
calculations  in  Fig.  9.  The  average  current  density  of  the  reference 
case  is  688  mAcnrr2.  For  the  three  cases  shown  in  Fig.  9,  the  aver¬ 
age  current  density  is  calculated  as,  760,  726,  and  563  mAcnrr2  for 
kc/kcjtd  =  1-0,  0.4  and  0.1,  respectively.  We  observe  approximately 
a  10%  improvement  for  kclkc_std  =  1.0.  Even  for  the  case  where 
kclkc_std  =  0.4,  about  a  6%  improvement  is  achieved.  It  is  worth  not¬ 
ing  that  the  improvement  in  the  average  current  density  is  achieved 
for  the  same  amount  of  cathode/electrolyte  material,  by  introduc¬ 
ing  a  non-uniform  cathode-electrolyte  interface  structure.  In  the 
case  where  /<c/kc_std  =  0.1 ,  however,  the  current  density  is  decreased 
by  18%. 

Fig.  10  shows  the  distributions  of  the  ion  potential,  activa¬ 
tion  overpotential,  concentration  overpotential  and  charge  transfer 
current  in  the  domain  of  the  optimization  calculation  for  the 
case  where  kclkcstd  =  1.0.  The  electron  potential  in  the  cathode  is 
assumed  uniform  in  this  study.  Therefore  the  potential  difference 
between  the  ion  and  the  electron  is  large  in  the  region  where  the  ion 
potential  is  high.  In  Fig.  10(a),  the  ion  potential  is  high  at  the  bot¬ 
tom  of  the  groove.  The  charge  transfer  current,  Itpb,  in  Fig.  10(d) 
takes  large  value  in  this  region  accompanied  by  high  activation 
loss  shown  in  Fig.  10(b).  As  a  result  of  intensive  oxygen  consump¬ 
tion,  the  concentration  overpotential  becomes  relatively  high  in 
the  same  region.  On  the  other  hand  at  the  convex  side  of  the  elec¬ 
trolyte,  ion  potential  is  relatively  low  because  of  the  ion  transport 
loss  in  the  dense  electrolyte.  Therefore  the  driving  force  of  the 
reaction  is  relatively  weak  compared  to  the  groove  bottom.  Never¬ 
theless  moderate  charge  transfer  current  is  observed  in  the  vicinity 
of  the  interface  even  at  the  top  region  of  the  convex  part  where  the 
concentration  overpotential  is  negligible.  Consequently  the  elec- 
trochemically  active  region  spreads  near  the  cathode-electrolyte 
interface. 

4.2.2.  Effect  of  complexity  parameter 

The  effects  of  the  complexity  parameter,  y  in  Eq.  (21),  upon 
the  optimized  structure  are  examined.  Optimization  calculations 
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are  conducted  for  three  different  values  of  the  complexity  param¬ 
eter  and  the  results  are  shown  in  Fig.  11.  Calculation  conditions 
other  than  the  complexity  parameter  are  the  same  for  the  three 
cases.  It  is  confirmed  that  different  optimization  results  can  be 
obtained  by  tuning  the  value  of  this  parameter.  In  general,  smaller 
values  of  the  complexity  parameter  are  seen  to  yield  optimal 
configurations  with  more  “complex”  geometry  [21].  The  shapes 
presented  in  Fig.  1 1  seem  to  follow  this  trend,  showing  a  relatively 
flat  cathode-electrolyte  interface  for  the  largest  value,  y  =  0.05. 
The  average  current  density  is  707,  726  and  694  mA  cm-2  for 
y  =  0.005,  0.02  and  0.05,  respectively.  The  results  indicate  that  cell 
performance  improvement  is  not  a  monotonous  function  of  the 
complexity  parameter. 

5.  Conclusions 

To  explore  potential  for  the  power  density  enhancement  of 
solid  oxide  fuel  cells  by  controlling  the  electrode-electrolyte  inter¬ 
face,  we  performed  two-dimensional  numerical  simulations  of  a 
SOFC,  where  a  level  set-based  topology  optimization  technique 
was  applied  for  the  local  optimization  of  the  interface  shape.  The 
following  conclusions  were  obtained. 

(1 )  The  level  set-based  optimization  method  was  successfully  cou¬ 
pled  with  the  SOFC  simulation  based  on  a  microscale  model 
for  the  first  time.  The  optimization  method  was  applied  to 
the  cathode  side  of  a  cell  to  find  the  optimized  shape  of  the 
cathode-electrolyte  interface  that  maximizes  the  current  den¬ 
sity  under  a  given  terminal  voltage.  Converged  solutions  were 
obtained  after  sufficient  iterations.  The  solutions  were  con¬ 
firmed  to  be  independent  of  the  initial  conditions. 

(2)  The  numerical  results  show  that  the  cell  performance 
can  be  improved  by  applying  non-flat  design  to  the 
cathode-electrolyte  interface  for  the  same  amount  of  cath¬ 
ode/electrolyte  materials.  In  a  standard  case  set  in  this  study, 
kclkc_std  =  \.0  and  y  =  0.02,  a  wavy  interface  is  obtained  as  an 
optimized  shape.  Approximately  a  10%  improvement  of  the 
current  density  is  observed  for  this  case  compared  to  a  con¬ 
ventional  flat  interface.  The  shape  of  the  interface  can  be 
tuned  by  changing  the  complexity  factor  of  the  optimization 
calculation. 

(3)  The  numerical  results  show  that  the  optimized  shape  of  the 
cathode-electrolyte  interface  varies  depending  on  the  simu¬ 
lation  conditions.  When  the  gas  diffusion  resistance  in  the 
cathode  is  large,  the  optimized  interface  tends  to  become  flat. 
An  optimized  interface  is  determined  as  a  result  of  combined 
effects  of  various  resistances  closely  related  to  the  microstruc¬ 
ture  of  the  electrode. 
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Appendix  A.  Preliminary  experiment  with 
cathode-electrolyte  interface  modification 

The  results  of  the  optimization  simulation  show  the  possibil¬ 
ity  of  the  cell  performance  enhancement  by  cathode-electrolyte 


interface  modification  in  mesoscale.  A  preliminary  experiment  is 
conducted  to  see  the  fundamental  effects  of  the  interface  modifi¬ 
cation. 

Fig.  Al  shows  an  overview  of  the  experimental  apparatus.  A 
button  type  cell  is  used.  The  entire  test  section  is  placed  in  an  elec¬ 
tric  furnace  to  control  the  operation  temperature.  Fig.  A2  shows 
the  electrolyte  of  a  test  cell,  a  small  disk-shaped  YSZ  ceramic.  The 
diameter  and  thickness  of  the  electrolyte  are  20  mm  and  0.5  mm, 
respectively.  The  cathode  side  surface  is  partially  grooved  by  apply¬ 
ing  sandblasting  method  on  the  initial  flat  surface.  Two  types  of 
groove  shape  are  prepared,  and  the  details  of  their  geometries  are 
summarized  in  Table  Al.  The  area  enlargement  factor  is  defined  as 
the  ratio  of  the  electrode-electrolyte  interface  area  of  a  grooved 


Fig.  A2.  Top  view  of  cathode  side  electrolyte  with  grooves  (large  grooves). 
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a  Large  grooves  (LG)  b  Small  grooves  (SG) 


Fig.  A4.  Cross-section  view  of  grooved  portions. 


Table  Al 

Geometric  details  of  the  grooves. 


Depth  (|Jim) 

Width  (|Jim) 

Interval  (|jim) 

Area  enlargement  factor 

Large  grooves  (LG) 

210 

210 

220 

1.58 

Small  grooves  (SG) 

60 

60 

50 

1.72 

portion  to  a  corresponding  flat  portion.  The  electrolyte  surface  is 
flat  on  the  anode  side. 

Electrode  materials  are  prepared  following  the  procedure  noted 
in  Ref.  [20].  The  anode  material,  NiO-YSZ,  is  mixed  with  polyethy¬ 
lene  glycol,  screen  printed  on  the  electrolyte  and  sintered  at  1400  °C 
for  5  h.  The  cathode  material,  LSM,  is  also  mixed  with  polyethylene 
glycol  to  form  slurry,  which  is  screen  printed  on  the  electrolyte 
and  sintered  at  1150°C  for  5h.  As  can  be  seen  in  Fig.  A3(b),  the 
cathode  is  segmented  into  two  parts.  One  side  covers  the  flat  sur¬ 
face  and  the  other  covers  the  grooved  surface.  LSM  cathode  is  used 
because  the  effect  of  the  grooves  is  expected  to  be  prominent  when 
the  electrochemical  reaction  takes  place  only  in  the  vicinity  of  the 
interface.  Considering  that  the  ionic  conductivity  of  LSM  is  small, 
LSM  cathode  is  expected  to  maximize  the  difference  between  the 
two-segmented  cathodes.  Fig.  A4  shows  cross-sectional  views  of 
the  grooves  for  the  two  types  of  cell.  Power  generation  tests  using 
each  part  of  the  cathode  are  conducted  individually  to  obtain  the 
i-V  characteristics  and  results  are  compared.  The  current  density,  i, 
is  defined  as  the  measured  current  divided  by  the  projected  cathode 
area.  The  projected  area  is  evaluated  through  image  processing.  The 


experiments  were  conducted  at  various  operating  temperatures 
(800,850  and  900  °C). 

Fig.  A5  compares  the  i-V  curves  obtained  from  two  cathodes 
with  and  without  the  grooves  at  different  temperature.  Dry  hydro¬ 
gen  is  supplied  to  the  anode  side  as  fuel  while  dry  air  is  supplied 
to  the  cathode  side.  The  current  density  obtained  from  the  cath¬ 
ode  with  grooves  is  always  larger  than  that  from  the  flat  cathode 
in  both  groove  sizes  (LG  or  SG)  when  it  is  compared  at  the  same 
terminal  voltage.  Fig.  A5  shows  an  experimental  example  where 
the  cell  performance  is  enhanced  by  the  cathode-electrolyte  inter¬ 
face  modification.  If  the  current  density  at  the  terminal  voltage 
of  0.6  V  is  compared,  the  cathode  with  grooves  achieves  approx¬ 
imately  30-40%  larger  value  than  the  flat  cathode.  This  is  less  than 
the  area  enlargement  factor  in  Table  Al .  As  discussed  with  the  sim¬ 
ulation  results,  the  cell  performance  is  affected  by  many  factors, 
local  ion  potential,  local  electron  potential  and  gas  composition, 
etc.  The  degree  of  enhancement  by  the  interface  modification  is 
determined  as  a  result  of  the  combination  effects  of  those  factors. 
To  clarify  each  effect  in  experiments  is  an  important  work  left  for 
future  investigation. 
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(a)  Large  grooves  (LG) 
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(b)  Small  grooves  (SG) 
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Fig.  A5.  i-V  characteristics  of  the  cells  with  segmented  cathodes.  In  each  cell,  one  cathode  is  fabricated  on  the  grooved  electrolyte  with  the  other  cathode  on  the  flat  electrolyte. 
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